## AQI Analysis - Response to  Reviewer 3 Interpretation Requests

rm(list = ls())

library(xtable)

# rid of E
options(scipen=999)
### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")
load("ppp_cleaned")

sd(p2$aqio) # 63.83

# median of moderate is 75, median of Unhealthy is 175, median of Hazardous is 400 

# cacl multipliers in terms of percentage increase, 
s1 <- 100/sd(p2$aqio) # 1.56 sd
s2 <- 325/sd(p2$aqio) # 5.09 sd

########## ###########
## Now rerun the main models to extract aqi coefficients
########## ###########

# aqi and government satisfaction at the local level
l_satc <- lm(l_sat
             ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
             + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov1 <-cluster.vcov(l_satc, p2$day)  
coeftest(l_satc, m.vcov1) # coeftest (this will give you results using clustered standard errors)
se11 <- coeftest(l_satc, m.vcov1)

NWvcov <- NeweyWest(l_satc, lag = NULL, prewhite = FALSE)
sen11 <- coeftest(l_satc, NWvcov) 
# aqi and c_sat
c_satc <- lm(c_sat
             ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
             + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov25 <-cluster.vcov(c_satc, p2$day)  
coeftest(c_satc, m.vcov25) 
se22 <- coeftest(c_satc, m.vcov25)

NWvcov <- NeweyWest(c_satc, lag = NULL, prewhite = FALSE)
sen22 <- coeftest(c_satc, NWvcov) 

# aqi and the demand for more watchdog at the local level
l_watchc <- lm(l_watch
               ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
               + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov13 <-cluster.vcov(l_watchc, p2$day)  
coeftest(l_watchc, m.vcov13) # coeftest

se33 <- coeftest(l_watchc, m.vcov13)
NWvcov <- NeweyWest(l_watchc, lag = NULL, prewhite = FALSE)
sen33 <- coeftest(l_watchc, NWvcov) 

# aqi and the demand for more watchdog at the central level 
c_watchc <- lm(c_watch
               ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
               + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov14 <-cluster.vcov(c_watchc, p2$day)  
coeftest(c_watchc, m.vcov14) 

se44 <- coeftest(c_watchc, m.vcov14)
NWvcov <- NeweyWest(c_watchc, lag = NULL, prewhite = FALSE)
sen44 <- coeftest(c_watchc, NWvcov) 


# aqi and anti-west sentiment
anti_westc <- lm(anti_west
                 ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
                 + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov25 <-cluster.vcov(anti_westc, p2$day)  
coeftest(anti_westc, m.vcov25) 
se66 <- coeftest(anti_westc, m.vcov25) 

NWvcov <- NeweyWest(anti_westc, lag = NULL, prewhite = FALSE)
sen66 <- coeftest(anti_westc, NWvcov) 

# day pollute
day_pollutec <- lm(day_pollute
                   ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
                   + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov25 <-cluster.vcov(day_pollutec, p2$day)  
coeftest(day_pollutec, m.vcov25) 
se77 <- coeftest(day_pollutec, m.vcov25) 
 
NWvcov <- NeweyWest(day_pollutec, lag = NULL, prewhite = FALSE)
sen77 <- coeftest(day_pollutec, NWvcov) 

# aqi and life stisfaction
life_satc <- lm(life_sat
                ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ 
                + kids + married + income + gender + ccp + parade + national + age, data=p2)
m.vcov25 <-cluster.vcov(life_satc, p2$day)  
coeftest(life_satc, m.vcov25) 

se88 <- coeftest(life_satc, m.vcov25) 
NWvcov <- NeweyWest(life_satc, lag = NULL, prewhite = FALSE)
sen88 <- coeftest(life_satc, NWvcov) 

############ ###########
# Now multiple multipliers by coefficents from main models
############ ###########


a <- list(l_satc, c_satc, l_watchc, c_watchc, anti_westc)

s1.effects <- unlist(lapply(a, function(x) x$coefficients["aqi"]))*s1
s2.effects <- unlist(lapply(a, function(x) x$coefficients["aqi"]))*s2

s1.effects

s.effects <- cbind(s1.effects,s2.effects)
s.effectsT <- t(s.effects)

colnames(s.effectsT) <- c("Loc. Govt.", "Cent. Govt.", "Loc. Oversight", "Cent. Oversight",  "Anti-West")

rownames(s.effectsT) <- c("Moderate -> Unhealthy", "Moderate -> Hazardous")

s.effectsT <- round(s.effectsT,digits = 3)


## TABLE A21
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/category_move_sims.tex")
print(xtable((s.effectsT)),floating=FALSE,font.size="tiny",latex.environments="center")
sink()
